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The Kuramoto model describes a system of globally coupled phase-only oscillators with different 
natural frequencies. The model in the steady state exhibits a phase transition as a function of 
the coupling parameter between a low-coupling incoherent phase in which the oscillators oscillate 
independently and a high-coupling synchronized phase in which they oscillate with a common fre- 
quency. Here, we consider a uniform distribution for the natural frequencies, for which the phase 
transition is of first oder. We study how the system close to the phase transition relaxes in time 
to the steady state while starting from an initial homogeneous state. In this case, the relaxation 
process occurs as a step-like jump in the order parameter, and has been previously interpreted as 
suggesting the existence of metastable states. By analyzing (i) the distribution of the jump time, 
and (ii) the distribution of the order parameter in the steady state, we provide evidence to rule out 
the existence of any metastable states. 

PACS numbers: 05.45.Xt 
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Systems of coupled oscillators of different natural fre- 
quencies often exhibit collective synchronization in which 
a finite fraction of the oscillators oscillates with a com- 
mon frequency. Examples include groups of fireflies flash- 
ing in unison [l|, Q , networks of pacemaker cells in the 
heart 0, 0] , superconducting Joscphson junctions 0, @ , 
and many others. Understanding the nature and emer- 
gence of synchronization from the underlying dynamics 
of such systems is an issue of great interest. A paradig- 
matic model in this area is the so-called Kuramoto model 
involving globally-coupled oscillators Q • Although stud- 
ied extensively in the past, the model continues to raise 
new questions, and has been a subject of active research; 
for reviews, see (HQ. 

The model consists of N phase-only oscillators labeled 
by the index i = 1,2, ... ,N. Each oscillator has its 
own natural frequency Wj distributed according to a given 
probability density g(w), and is coupled to all the other 
oscillators. The phase of the oscillators evolves in time 
according to Q 
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where 6i, the phase of the zth oscillator, is a periodic 
variable of period 2ir, and K > is the coupling constant. 

The Kuramoto model has been mostly studied for a 
unimodal g(uj), i.e., one which is symmetric about the 
mean frequency u> = £1, and which decreases monoton- 
ically and continuously to zero with increasing \u> — Q\ 
[H, Then, it is known that in the limit N — > oo, the 
system of oscillators in the steady state undergoes a con- 
tinuous transition at the critical threshold K c = 2/7rg(0). 
For K < K c , each oscillator tends to oscillate indepen- 
dently with its own natural frequency. On the other 
hand, when K > K c , the coupling synchronizes the 



phases of the oscillators, and in the limit K oo, they 
all oscillate with the mean frequency fl. 

The degree of synchronization in the system at time t 
is measured by the complex order parameter 



1 



r(t) = r(t)e*W = - 
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with magnitude r(t) and phase ip(t), in terms of which 
the time evolution (TJJ may be written as 



dt 



Kr(t) sin(V>(£) - 
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Here r(t) with < r(t) < 1 measures the phase coher- 
ence of the oscillators, while ip(t) gives the average phase. 
When K is smaller than K c , it is known that r(t) while 
starting from any initial value relaxes in the long-time 
limit to zero, corresponding to an incoherent phase in 
the steady state. For K > K c , on the other hand, r(t) 
grows in time to asymptotically saturate to a non-zero 
value r(K) < 1 that increases continuously with K. 

For a unimodal g(uj), the relaxation of r(t) to steady- 
state is exponentially fast for K > K c . For K < K c , 
however, the nature of relaxation depends on g(ui). When 
g(uj) has a compact support, r(t) while starting from any 
initial value decays to zero more slowly than any expo- 
nential as t — > oo [To|. When g(u) is supported on the 
whole real line, r(t) as a function of time is known only 
in particular cases. For example, for a Lorentzian g(oj), 
and a s pec ific initial condition, r(t) decays exponentially 
to zero For other choices of g(oj) in this class and for 
other initial conditions, the dependence of r(t) on time is 
not known analytically, and it has bee n sp eculated that 
r(t) is a sum of decaying exponentials [lfj. 

In the limit N — > oo, the state of the oscillator sys- 
tem at time t is described by the probability distribu- 
tion f(6,t,u>) that gives for each natural frequency w 
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the fraction of oscillators with phases between 6 and 
9 + dd at time t. The time evolution of f(6,t,ui) sat- 
isfies the continuity equation for the conservation of the 
number of oscillators with natural frequency lo, and is 
given by a non-linear partial integro-differential equation 
|8j. Recent analytical studies for a unimodal g(ui) (specif- 
ically, a Lorentzian) and for a bimodal g(ui) (sum of two 
Lorentzians) demonstrated by considering a restricted 
class of f(6, t, u>) that the time evolution in terms of the 
integro-differential equation may be exactly reduced to 
that of a small number of ordinary differential equations 
(ODEs) [HI, [3. Interestingly, this system of ODEs for 
the reduced system contains the whole spectrum of dy- 
namical behavior of the full system. However, the valid- 
ity of such a reduction for a general g(u>) is not known. 

A uniform <?(w) with a compact support does not qual- 
ify as a unimodal distribution. In this case, it is known 
that in the limit N — ¥ 00, the Kuramoto model in the 
steady state exhibits a first-order phase transition be- 
tween an incoherent and a synchronized phase at the 
critical coupling K c = 2/7rg(0) For finite N, study 
of the relaxation of an initial state with uniformly dis- 
tributed phases has demonstrated that for values of K 
around K c , the process occurs as a step-like jump in r(t) 
from its initial to the steady state value. This behav- 
ior has been interpreted as sug gest ing the existence of 
metastable states in the system [lj, [l5| . 

In this work, we investigate in detail the implication 
of the existence of the step-like relaxation. By perform- 
ing extensive numerical simulations, we analyze in detail 
the distribution P(t) of the jump time r, and also the 
distribution P(r) of the steady state value r of the or- 
der parameter. Our analysis suggests the absence of any 
metastable states in the system. 

The numerical simulations involved integrating Eq. ([3]) 
by a 4th-ordcr Runge-Kutta algorithm. In simulations, 
we monitored the time evolution of r(t) while starting 
from an initial state with O^s independently and uni- 
formly distributed in [— 7r,7r] and w^'s independently and 
uniformly distributed in [—2,2] (so that g(u>) = 1/4). 
Thus, initially, r(t) w 0. We call such an initial state a 
homogeneous state. 

For K = 2.62, close to the critical value K c = 
2/715(0) w 2.55, Fig. QJa) shows the temporal evolu- 
tion of r{t) for 8 different realizations of the initial ho- 
mogeneous state for N = 1000, showing that r(t) jumps 
in a step at times ~ 100 from its initial to its steady 
state value ss 0.82 (with fluctuations). For the realiza- 
tion shown in Fig. D^b), however, the jump to the steady 
state value takes place at time sa 2000, while for the re- 
alization in Fig. [He), r(t) stays close to the initial value, 
and the jump does not take place in the time duration 
shown. The average of r(t) over the realizations in (a)-(c) 
is shown in Fig. [ljd), which naturally exhibits a step-like 
relaxation. Similar behavior of the average r(t), where 
the averaging is over a few realizations, was reported in 
fl4| . Figure QJe) shows the time evolution of r(t) when 
averaged over a large (~ 1000) number of realizations, 
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FIG. 1: (a) Evolution of r{t) in time for 8 different realizations 
of the initial homogeneous state. Here, K = 2.62, while the 
critical coupling is K c ~ 2.55. The plot shows that for these 
realizations, r(i) jumps in a step at times ~ 100 from its 
initial to its steady state value, (b) and (c) also show the 
time evolution of r{t) for two other realizations of the initial 
state: for the realization in (b), r(t) shows a jump to the 
steady state value at time « 2000, while for the one in (c), 
r(t) remains close to its initial value and does not relax to 
the steady state value in the time shown, (d) and (e) show 
respectively the time evolution of r(t) when averaged over the 
realizations in (a)-(c), and over typical 1000 realizations. The 
data are obtained from numerical integration of the iV-body 
dynamics © with N = 1000. 



which however does not show any step-like relaxation. 

Figure Q] illustrates that the relaxation of r(t) occurs 
as a step-like jump, and that there is a wide distribution 
of the jump time r. Thus, it is worthwhile to study the 
distribution P{t). Our numerical data, displayed in Fig. 
[2j show that P(t) has a power-law tail; the scaling col- 
lapse of the data, shown in Fig. [3l suggests the following 
scaling form: 



P{t) ~ v Nf 



(4) 



where f(x) is the scaling function. 

In the power-law regime of P{t), let us write 

P{t) OCT"", 



(5) 



where a is the decay exponent. Then, in this long-time 
regime, r(£) at time t will have a value close to one pro- 
vided that the jump from the initial close-to-zero value 
has already taken place by this time; otherwise, r(t) w 0. 
It then follows that the average order parameter, (r)(t), 
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FIG. 2: Distribution P(r) of the time r at which r(t) jumps 
from a close-to-zero value corresponding to an initial homoge- 
neous state to the final steady state value (see Fig. QJb) for an 
example of such a jump). The data are obtained from numer- 
ical integration of the JV-body dynamics (0. Here, K — 2.62. 
The inset shows power-law fits for the tail of P(t) according 
to Eq. (J3J), giving values 2.2, 2.7, and 2.8 for system sizes 
N = 1000, 2000 and 4000, respectively. 
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FIG. 3: Scaling collapse of the data for the distribution P(t) 
shown in Fig. [2] 



in this regime will be 

(r)(t) ~ / dr P(r) ~ 1 -t 1 -". (6) 
Jo 

Equivalently, one has 

l-{r)(t)~t-P, (7) 

with 

(3 = a-l. (8) 




FIG. 4: (a) Power-law fit to the tail of P(r), according to Eq. 
([5]). (b) Power-law fit to the quantity 1 — (r)(t), according to 
Eq. Q . The data are obtained from numerical integration of 
the iV-body dynamics © with N = 1000, and for K = 2.62. 
The average is computed over ~ 1000 realizations. 



This prediction is easily checked from our simulation re- 
sults. Figure SJb) shows the power-law fit to the tail 
of P(r) for N = 1000, giving a sa 2.2. On the other 
hand, Fig. S[b) shows the power-law fit to the quantity 
1 — (r)(t) at long times, giving (3 ps 0.9 in Eq. ([7]), while 
the value predicted by Eq. is 1.2. One may hope to 
obtain a closer agreement with better statistics. 

To probe the existence of metastable states in the 
system, we obtained the probability distribution of the 
steady state value r of the order parameter. Figure [5] 
shows that with increasing system size, the distribution 
is peaked at a value close to one, while the probability 
to obtain any other value is negligibly small. We have 
checked that the distribution P(r) shown in the figure 
does not evolve in time, and thus, corresponds to the 
steady state. The form of P(r) confirms the absence 
of metastable states in the thermodynamic limit, whose 
presence would have otherwise led to peaks in P(r) at 
values close to zero. 

Based on the above analysis, we thus conclude that 

1. For uniformly distributed frequencies and a homo- 
geneous initial state, the order parameter in the 
Kuramoto model shows relaxation as a step-like 
jump from its initial value (« 0) to its steady state 
value (~ 1) for values of K close to the phase tran- 
sition (K>K C ). 

2. There is a wide distribution of the time r of occur- 
rence of the jump, and its distribution P(t) has a 
power-law tail. 

3. Hence, averaging r(t) over a few realizations nat- 
urally shows a step-like relaxation, and may erro- 
neously lead to the conclusion that the system dur- 
ing the relaxation gets stuck in metastable states. 

4. In order to capture the whole range of values of r 
as suggested by its fat-tailed distribution, it is thus 
important that r(t) is averaged over a very large 
number of realizations. We demonstrated that in 
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doing so, r(t) does not show any step-like relax- 
ation, but instead a smooth evolution towards the 
steady state. This is consistent with the absence 
of any metastable states in the system, and is fur- 
ther corroborated by our results on the distribution 
P(r) of the steady state order parameter, which is 
sharply peaked around a value close to 1 in the 
thermodynamic limit. 

We obtained a scaling form for the distribution P(t). 
We also predicted and verified numerically a power-law 
growth of the averaged order parameter at late times. 
However, we do not have an understanding of the relax- 
ation mechanism, and an investigation in that direction 
is in progress. 



FIG. 5: Starting from an initial homogeneous state, the figure 
shows the distribution P(r) of the steady state value r of 
the order parameter. The data are obtained from numerical 
integration of the iV-body dynamics ((3| for the values of N 
shown in the figure, and for K = 2.62. 
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